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ABSTRACT 

In previous works we proposed the Reverse Zeldovich Approximation (RZA) method, 
which can be used to estimate the cosmological initial conditions underlying the galaxy distri- 
bution in the Local Universe using peculiar velocity data. In this paper, we apply the technique 
to run constrained cosmological simulations from the RZA-reconstructed initial conditions, 
designed to reproduce the large-scale structure of the Local Universe. We test the method with 
mock peculiar velocity catalogues extracted from a reference simulation. We first reconstruct 
the initial conditions of this reference simulation using the mock data, and then run the recon- 
structed initial conditions forward in time until z = 0. We compare the resulting constrained 
simulations with the original simulation at z = to test the accuracy of this method. We also 
compare them with constrained simulations run from the mock data without the addition of 
RZA, i.e. using only the currently established constrained realizations (CR) method. Our re- 
simulations are able to correctly recover the evolution of the large-scale structure underlying 
the data. The results show that the addition of RZA to the CR method significantly improves 
both the reconstruction of the initial conditions and the accuracy of the obtained constrained 
resimulations. Haloes from the original simulation are recovered in the re- simulations with an 
average accuracy of « 2 Mpc/h on their position and a factor of 2 in mass, down to haloes 
with a mass of « 10 14 M o //z. In comparison, without RZA the re- simulations recover only the 
most massive haloes with masses of « 5 • 10 14 M Q /h and higher, and with a systematic shift 
on their position of about ^10 Mpc/h due to the cosmic displacement field. We show that 
with the additional Lagrangian reconstruction step introduced by the RZA, this shift can be 
removed. 

Key words: cosmology: theory - dark matter - large-scale structure of Universe - galaxies: 
haloes - methods: numerical 



1 INTRODUCTION 

While numerical simulations have developed into a cornerstone of 
studying the large-scale structure (LSS) of the Universe, there is 
still a long way to go towards reconciling the predictions drawn 
from such cosmological simulations with observational data. On 
the one hand, the general properties of large-scale matter clustering 
and interacting are now very well understood. This process is typ- 
ically simulated by generating a random realization of the primor- 
dial density fluctuations of the Universe, and then integrating it for- 
ward in time using Af-body techniques and the physics defined by 
the concordance ACDM model. These simulations produce model 
universes which are statistically in good agreement with the ob- 
served LSS. On the other hand, the region best studied observa- 
tionally - the Local Universe - shows many properties that cannot 
be directly modelled with such random realizations. In order to ex- 
plain the properties and dynamics of the Local Group, i.e. the Milky 



Way and Andromeda galaxies and their satellites, one has to study 
its formation history, which seems to be tightly connected to the pe- 
culiar alignment of the large-scale structure in the Local Universe, 
such as the Local Void, the Local Supercluster (LSC), and farther 
away structures like the Great Attractor (GA) and Perseus-Pisces 
cluster. This also includes studies of the Local Universe's velocity 
field, since it directly traces the gravitational potential and therefore 
the total matter distribution. 

The ansatz of the CLUES projecfl, which provides the frame- 
work for the study presented here, is to conduct constrained 
simulations that are able to reprod u ce the LSS of the ob served 
Local Universe (Klypinetal. 2003; Gottlo ber et all ko 10). Such 
simulations serve as an ideal laboratory for studyin g structure 
formation in our cosmological neighbourhood dLibeskind et al.l 

1 www.clues-project.org 
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l20ld : iForero-Romero et al"1l201 ll: iKnebe et al.ll201 lal fbh. They are 
constructed by constraining the initia l conditions using observa- 
tional data faoffman & Ribald Il99ll : iBistolas & Hoffman! Il998l : 
IZaroubi et alj Il995l Il999h . in particular we will use the mea- 
sured radial pecul i ar velocities of galaxies in the Local Universe 
(iTullv etalJ[2008l . 120091 : ICourtois et alll2012h as the input data. 
To construct initial conditions for t hese simulations, the c on- 
strained realizations (CR) algorithm dHoffman & Ribald 1 1 99 ll) is 
used, which combines a Bayesian reconstruction from the data with 
a random component created with a conventional initial conditions 
generator. 

From previous CLUES simulations (Klv pin etalJ [2003; 
iGottlober et aljEoiob we know that the technique of using pecu- 
liar velocities as constraints has its limits. These sim ulations used 
radial velocities from the now out dated MARK III Jwiliick et alJ 
Il997h . SBF (iTonrv et alfeOOlh . and lKarachentsev et all J2004h cat- 
alogues. To obtain reasonable reproductions of the Local Uni- 
verse's structure, these velocity constraints needed to be comple- 
mented by additional densi ty constraints, which were d rawn from 
X-ray selected cluster data (IReiprich & Bohringerl l2002). But even 
in this case only the few most massive clusters of the Local Uni- 
verse, namely the LSC, the GA, and for larger boxsizes the Coma 
and Perseus-Pisces superclusters, were robust features appearing 
in the constrained realizations, with smaller scales essentially un- 
constrained and dominated by the random component. It is im- 
portant to understand whether these limitations stem from insuf- 
ficient accuracy of the observational data used as constraints, the 
method employed to construct the constrained initial conditions, 
or some fundamental physical limitation. Through a coll aboration 
of CLU E S with the observational Cosmicflows program (Courtoisl 
l2011al lbl: ICourtois & TuUvl 120121 : ITullv & Courtoisl 1201 2^ it "will 
become possible to use significantly higher-quality data with many 
more peculiar velocity datapoints out to higher distances, with bet- 
ter sky coverage and much smaller observational errors, for gen- 
erating constrained simulations. Here, we want to investigate how 
much increase in accuracy we can expect from constrained simula- 
tions set up with this new data compared to what we have now, and 
how the method itself of setting up constrained initial conditions 
can be improved in order to optimally utilize the additional infor- 
mation contained in the data. With these studies, we hope to pave 
the way for a new generation of accurate constrained simulations 
that will be produced by our collaboration during the next years, to 
provide a powerful framework in which the Local Universe can be 
studied in detail. 

This work is the third in a series of papers on this subject. 
In the first paper (Doumler et al. 2012a, from here on Paper I) we 
presented the Reverse Zeldovich Approximation (RZA) method, 
which significantly increases the quality of reconstructed initial 
conditions obtained from peculiar velocity data. This is accom- 
plished by a Lagrangian reconstruction of the primordial density 
distribution underlying the observed velocity field, essentially shift- 
ing the data back in time and thus provide better constraints for 
initial conditions than the original dataset observed at z = 0. In the 
second paper (Doumler et al. 2012b, from here on Paper II) we then 
investigated the impact of observational errors on the RZA method. 
Here, we study how well constrained simulations can reproduce the 
underlying universe, if their initial conditions were constructed by 
employing the RZA method. For this we set up a detailed test by 
using realistic mock peculiar velocity data drawn from a test sim- 
ulation snapshot at z = 0. We use the data to generate constrained 
initial conditions at some early initial redshift Zinit and run them for- 
ward again with an A/-body code. This evolved re-simulation is then 



compared to the original simulation at z = 0. We do this for both 
the previous CLUES method and our new RZA method to compare 
by how much the simulation accuracy improves by performing La- 
grangian reconstruction on the data. We also compare two mocks 
of different quality, to estimate how much is gained by using more 
accurate data. 

The outline of this paper is as follows. In Section[2] we briefly 
review our method of generating constrained initial conditions from 
the data, describe the setup of this test, and present the set of re- 
simulations we conducted. In Section [3] we study the accuracy of 
these re- simulations compared to the original reference simulation 
and present our findings. We summarize and discuss our results in 
Section|4] 



2 METHOD AND TEST SETUP 
2.1 Initial conditions with RZA 

Our method of RZA reconstruction and subsequent generation of 
constrained initial conditions is described in detail in Paper I; we 
give only a brief summary here. 

We start with a set of datapoints, i.e. radial peculiar veloci- 
ties v r at discrete positions r at z = 0. We first apply a grouping 
procedure to the data in order to "linearize" it, i.e. to remove virial 
motions and other small-scale interactions inside galaxy groups and 
clusters and to produce a data set that traces the coherent large-scale 
velocity field. We then reconstruct the three-dimensional peculiar 
velocities v(r) at positions r by using the Wiener Filter (WF). The 
WF produces an estimate u WF (r) based on the correlation function 
given by the assumed prior model, which is defined through the 
cosmological parameters and power spectrum P(k). It also filters 
out noise due to observational errors from the data. In order to con- 
struct cosmological initial conditions at some early redshift z ini t, we 
need to obtain a suitable set of constraints. We do this with the RZA 
method: we estimate the initial position x init at z ini t of our peculiar 
velocity field tracers, which are located at r at z = 0, with 



H f 



(1) 



The reconstructed initial positions x init are not to be understood as 
actual positions of the observed galaxies at early times - at typical 
initial redshifts z ini t, no galaxies have formed yet. They are rather 
interpreted as a way to trace the estimated initial peculiar velocity 
field at Zinit- In order to construct the full set of constrained initial 
conditions, we then shift the original datapoints v r "back in time" 
to Xi n it and use them as constraints to con struct a constrained real- 
ization with the lHoffman & Ribakl |l991) method. We only change 
the position of the constraints from r to x init , but preserve the am- 
plitude of the velocity and the direction of the component that is 
constrained (note that it is not the radial direction anymore w.r.t. 
the observer because the position has changed). The constrained 
realization is then obtained by evaluating 



S CR (r) = S RR (r) + (S(r) Ci ) (c iCj )- 1 ( Cj - cj) 



(2) 



where 6 RR is an independently generated random realization, c t are 
the RZA-shifted datapoints, c t are the corresponding values of the 
same quantities in the random realization, and the angled brackets 
denote the values of the correlation functions of the different quan- 
tities, defined by P(k). In order to construct an actual set of initial 
conditions for an Af-body simulation, we scale S RR to Zinit and solve 
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Figure 1. Slice through the BOX 160 constrained simulation from where 
the mock catalogues have been extracted. The cross is placed at the posi- 
tion tmw of the mock observer. The dashed circles illustrate the mock data 
volume for mocks with Rmax = 30 Mpc/h and 60 Mpc/h, respectively. 

for the displacement field, i/f(r) = -V _1 £(r). We then place parti- 
cles on a grid with displacements i/f according to the e stablished 
Zeldovich- approximation method (Efstat hiou et al.|[l985l) . 



2.2 Mock data 

As we already did in Papers I and II, we use the BOX 160 simulation 
conducted by the CLUES project as the source of our mock galaxy 
peculiar velocity data. Again, we refer the reader to these papers for 
details. The BOX160 is a constrained simulation of the Local Uni- 
verse with a boxsize of 160 Mpc/h, set up with the WMAP3 cos- 
mological parameters. The simulation contains a large-scale struc- 
ture resembl i ng the observed Local Universe (see Figure Q~] and 
ICuesta et alJ bOllh ). The simulation contains a configuration of 
objects corresponding to the LSC, GA, Coma, and Perseus-Pisces 
clusters (labelled in Figure [TJ, and a Local Group (LG) candidate. 
We select the position of this LG object (marked as a white cross 
in Figure [TJ as the mock observer and generate from there real- 
istic mock observational catalogues of galaxy peculiar velocities. 
A detailed description of how the mocks are built can be found in 
Paper I. Here we want to mention again that the mocks realisti- 
cally reproduce features of real observational catalogues, such as a 
limited distance, knowledge of only the radial component v r of v, 
sparse sampling, and observational errors due to inaccurate galaxy 
distance measurements. In this paper, we concentrate on the partic- 
ular mocks C30_10 and E60_10, which are designed to mimic the 
current Cosmicflows- 1 catalogue and the upcoming Cosmicflows- 
2 catalogue, respectively. After grouping, the catalogue C30_10 
contains 588 radial peculiar velocity datapoints within a relatively 
small radius of 30 Mpc/h from the mock observer. This is a simi- 
lar quality like the datasets that were previously used to construct 
constrained simulations. The E60_10, on the other hand, contains 



7632 datapoints within 60 Mpc/h, and models the quality of the 
upcoming new Cosmicflows-2 dataset, which we plan to use for 
constrained simulations of the Local Universe in future work. The 
data zone radii of 30 and 60 Mpc/h, respectively, are marked in 
Figure Q] with dashed white circles. 

Since we want to study specifically how well realizations of 
cosmological initial conditions can be constrained with peculiar ve- 
locity data, and how the RZA method performs in this context, in 
this work we do not use any other types of constraints such as clus- 
ter density constraints for our re- simulations. 

2.3 The set of constrained realizations 

For each of the two mocks, C30_10 and E60_10, we construct sev- 
eral constrained realizations of cosmological initial conditions. We 
use the method outlined in Section |2TI which combines the CR 
method with RZA reconstruction. In the following, we refer to this 
procedure as "Method II". Furthermore, we also generate initial 
condi tions with the method previously used for CLUES simula- 
tions dKlvpin et alJ2003l : lGottlober et al.l201oh . which we call here 
"Method I". It consists of using the peculiar velocity datapoints at 
Z = directly as constraints for Eq. (|2j, omitting the RZA shift. The 
main drawback of Method I is that it treats the peculiar velocities as 
though linear theory would be valid at all scales at z = 0, neglect- 
ing all higher-order effects such as the cosmological displacement 
field i/f. This leads not only to a poorer reconstruction quality, but 
also to a systematic position error of the clusters recovered in such 
simulations (compared to their observed counterparts). Typically, 
the object's positions will be off by the amplitude of i//, which is 
about 10 Mpc/h on average at z = 0. These shifts were observed in 
all previous CLUES simulations. We showed in Paper I that RZA 
can compensate for the shifts; here, we want to demonstrate how 
this improved method affects the outcome of evolved constrained 
simulations at z = 0. 

Having two different mock catalogues and two different meth- 
ods to generate ICs, we also want to study the impact of the random 
component 6 RR in Eq. (|2]). The peculiar velocity constraints c t af- 
fect only large scales from « 5 Mpc/h upwards, and only in regions 
of the box well covered by the data; all other structures emerging 
in the constrained simulation will have their origin in the particu- 
lar realization of the random component. We therefore expect that 
the random seed has a large impact on the outcome of the simula- 
tion. Varying the seed while keeping all other parameters such as 
the constraints constant allows us to estimate how robustly struc- 
tures that are constrained by the data are actually recovered in the 
constrained simulations. 

For each mock-method combination, we created six different 
realizations with different seeds for the random component. We 
have therefore a set of 24 different realizations of initial conditions. 
We use this set to test our method on scales smaller than the box. In 
all cosmological simulations - constrained as well as unconstrained 
- due to periodic boundary conditions the dynamics on the scale 
of the box is incorrect. Therefore, one must expect that also in fu- 
ture constrained simulations based on observational data the bulk 
flow on scales of the order of the simulation box will be incorrect. 
We construct the initial conditions on a regular cubic grid with a 
resolution of N = 25 6 3 and a boxsize of L = 160 Mpc/h, match- 
ing the boxsize of the "source simulation" BOX16C0. All realiza- 

2 The original BOX 160 simulation was performed with the ART code 
(Krav tsov et al.lll997h and a resolution of N = 1024 3 . However, since it 
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tions were constructed with our newly developed numerical code 
ICeCoRe (see Paper I for details). For the 24 different realizations, 
we assumed the following naming convention. We abbreviate the 
six different seeds as A through F. We then add the method and 
seed numbers to the end of the mock name, so that for example 
C30_10_II_A refers to the first out of six realizations that were con- 
structed with constraints from the C3CL10 mock using Method II. 
We then ran each of the generated constrained realizations of initial 
conditions forw ard until z = with the simulation code Gadget-2 
(Springel 2005), using collisionless particles only (no SPH parti- 
cles) with a resolution of N = 25 6 3 particles. We also used the 
same cosmological parameters and initial power spectrum P(k) that 
was used for the BOX 160 simulation, to be fully consistent on the 
assumed cosmological model. In fact, assuming a different cosmol- 
ogy would change the result. For example, increasing substantially 
the normalization of the power spectrum without changing the con- 
straints leads to a faster evolution. Instead of a local group like ob- 
ject one would find at the same place a massive large halo of about 
the total mass of the group. 

3 RESULTS 

3.1 Scatter and mass function 

Figure [2] shows a cell-to-cell comparison between the evolved con- 
strained resimulations obtained from the C30_10 mock and the 
original field of BOX 160 in the constrained volume (out to 30 
Mpc/h distance from the mock observer), both for Method I (with- 
out RZA) and Method II (with RZA). Figure [3] displays the same 
for the resimulations from the E60_10 mock. It can be seen that 
all mocks and methods produce density and velocity fields that 
are well-correlated with the original BOX 160 universe at z = 0. 
However, especially for the more accurate E60_10 mock, the added 
RZA reconstruction that was used to obtain the initial conditions in 
Method II significantly improves the correlation over the one ob- 
tained with Method I (without RZA). For the E60_10 mock, the 
peculiar velocity field is reproduced accurately down to an accu- 
racy of about 1/4 of its total standard deviation with Method II 
that uses RZA. This is an impressive improvement of the recon- 
struction quality. From the correlation of the density fields, it can 
be seen that Method I, which does not use RZA, significantly un- 
derestimates the total density of the constrained region, which is 
overdense compared to the cosmic mean. The RZA method, on the 
other hand, accurately reproduces this total overdensity, for both 
the C30_10 and E60_10 resimulations. The total overdensity as well 
affects the abundance of dark matter haloes. Figure |4] shows the 
binned dark matter halo mass functions of all realizations obtained 
from the C30_10 mock catalogue, both with and without RZA. 
The mass functions of the individual realizations are shown with 
coloured lines and points, their average is represented by the thick 
dashed black line, and the actual mass function of the BOX 160 in- 
side the 30 Mpc/h region is shown with the solid dashed black line. 
This is an overdense region in the BOX160 simulation: the mass 
function of this region lies significantly above the cosmic mean 
(grey line). The volume contains about 2.5 times more haloes with 
masses > \0 13 M Q /h than a region of mean cosmic density would 

is impractical to run a large number of re- simulations with such a resolu- 
tion, we limited ourselves to N = 256 3 and also re-ran the BOX 160 with 
Gadget-2 and this lower resolution to serve as the reference simulation. This 
way we can stay fully consistent on the parameters of all simulations. 
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Figure 2. Cell-to-cell comparison between the evolved constrained resimu- 
lations at z = for and the original BOX 160 simulations. Top row: velocity 
field inside the mock volume (all three components were concatenated); 
Bottom row: density field. All fields were smoothed with a 5 Mpc/h Gaus- 
sian. 

with the given cosmological model and boxsize. The constrained 
resimulations without RZA only partially follow that behaviour. 
Their average mass function is also above the cosmic mean, but still 
systematically underestimates the actual halo abundance in the cor- 
responding region of the original BOX 160 simulation. On the other 
hand, the resimulations where RZA was used for constructing the 
initial conditions follow the mass function of BOX 160 much more 
closely. 



3.2 The BOX160 universe 

The large-scale structure of the "local universe" for our test case, 
i.e. the region around the mock observer in the BOX 160 simu- 
lation, is shaped by several dominant structures that to some de- 
gree resemble the observed Local Universe. The "Milky Way" ana- 
logue halo that was chosen as the position of the mock observer 
(X = 75; Y = 64; Z = 80) lies in the vicinity of a "Virgo" clus- 
ter with a virial mass of 3.25 x 1O 14 M , which is embedded in a 
thick filament parallel to the X-axis. The local flow of the mock 
observer is directed towards this struc ture, resembling the observed 
Virgocentric infall dTullv et al.ll2008l) . On a larger scale, the whole 
region is dominated by a flow towards the "Great attractor" (GA), 
a massive structure at about X = 30, Y = 60 that lies outside of 
the 30 Mpc/h data zone of the C30_10 mock. This configuration 
can be seen in the top left map of Figure [5] which is a zoom-in 
from Figure Q] The BOX160 Virgo cluster is however not the most 
massive structure within 30 Mpc/h, as there is an even more mas- 
sive region at a distance of slightly less than 30 Mpc/h, lying in a 
direction in between Virgo and the GA, that contains two clusters 
with masses of 6.06 x 10 14 M Q /h and 5.20 x 10 14 MJh, respectively 
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Figure 4. Dark matter halo mass function inside the 30 Mpc//z mock volume for the six different realizations from C30.10 without RZA (top) and with 
RZA (bottom). The individual realizations are shown with coloured lines/points; their average mass function is the thick black dashed line. The original mass 
function of the BOX 160 in this volume is the thic k solid black lin e. The mass function for an average sub volume of radius 30 Mpc//z with the BOX 160 
parameters is shown in grey (theoretical model of Tin ker et al . 2008). 
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Figure 5. Density and velocity fields in a 10 Mpc//z thick slice at 87 < Z < 97 Mpc//z for the original BOX 160 simulation (top left) and three constrained 
resimulations without RZA (top right), with RZA (bottom left), and with RZA from the E60.10 mock using a larger data volume (bottom right). The white 
circle marks the constrained region for the C30.10 mock. The density is shown in logarithmic scale. The arrows are proportional to the amplitude of the 
peculiar velocity at each position. The prominent object inside the BOX160 data zone is the simulated Virgo cluster with a virial mass of 3.25 x 10 14 M Q /h. In 
all three constrained resimulations the same random seed was used. 



(their Z positions are just outside of the slice plotted in Figure [5). 
We associate them with the Centaurus and Hydra clusters. They in 
turn cause a significant infall flow towards them on the surround- 
ing structure. In the other direction (towards negative Y), there is 
also a large void that contributes to the shape of the large-scale 
velocity field by creating a push outwards of it, although this par- 
ticular feature in t he B OX 1 60 may not be as dominant as the ob- 
served Local Void dTuliv etal. 2008). Constrained resimulations of 
this test universe should be able to recover all these characteristic 
structures. BOX160 contains this specific configuration because as 
already mentioned it is itself a constrained simulation of the Local 
Universe. It is not entirely accurate, for example the virial mass of 
the BOX160 Virgo cluster is 2 - 3 times less than the estimated 



virial mass of its observed counterpart (iFouque et al.ll200lk but 
the described main characteristics of the large-scale structure are 
present. Ignoring the origin of BOX 160 for this test and taking its 
large-scale matter distribution as "given", we try to reproduce it in 
the constrained resimulations, which should retain the main char- 
acteristics that BOX 160 shares with the observed Local Universe 
in this "second pass" of the reconstruction-resimulation cycle. 

3.3 The re-simulated Virgo cluster 

The other panels of Figure|5]show three constrained realizations ob- 
tained from the C30_10 mock catalogue without RZA (top right), 
with RZA (bottom left), and from the larger-volume E60_10 mock 
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Figure 3. As Figure [2] but for constrained realizations from the E60.10 
mock with twice the data zone radius. 



catalogue with RZA, in the 10 Mpc/h thick slice that should contain 
the BOX 160 Virgo cluster. One can immediately notice that, while 
all three simulations faithfully reproduce the large-scale flow to- 
wards the GA, the resimulation that was constructed without RZA 
fails to create a Virgo cluster. Apparently, the mass of the BOX 160 
Virgo cluster of 3.25 x \0 14 M Q /h lies below the scale that could be 
reproduced with the non-RZA Method I of generating constrained 
initial conditions. Why then does the BOX160, which was created 
with the same method, have a Virgo cluster in the first place? There 
are three reasons: First, the observed Virgo cluster has a mass that 
is 2 - 3 times higher, which places it in the recoverable mass range. 
Second, we observe a scatter in the recovered masses and BOX160 
is a realisation with high Virgo mass. Third, in the BOX160, the 
Virgo cluster was additionally constrained by density constraints 
placed on its observed position. In our setup, we do not use such ad- 
ditional constraints. The non-RZA resimulations show some over- 
density in that region, and there is also the tendency of a flow in that 
direction, but there is not enough overdensity to form a cluster of 
comparable mass. On the other hand, both resimulations in Figure 
\5\ that include RZA reconstruction have a cluster at that position, 
which is similarly embedded in a thick filament, with a strong flow 
towards it from the position of the mock observer. 

In order to understand more clearly how robustly the Virgo 
cluster is recovered in the constrained resimulations we searched 
for it in all resimulations by searching for haloes that would be 
within 10 Mpc/h of the original BOX 160 Virgo's position and 
would have a mass of at least a 10 14 M Q /h. The result was that 
in all realizations using RZA, both from the C30_10 mock and 
the E60_10 mock, we could find a cluster corresponding to Virgo. 
These objects are listed in Table Q] On the other hand, we could not 
find such an object in any of the realizations created without RZA, 
which clearly displays that the RZA improves recovering the orig- 
inal structures in constrained simulations. The resimulated Virgo 



is not exactly at its BOX 160 position in the resimulations; the er- 
ror lies between 1.7 and 6.3 Mpc/h and varies with the random 
seed. The average position of all resimulated Virgos is only about 
2 Mpc/h away from its original position in BOX 160, so this fluctu- 
ation is probably due to the influence of the random modes rather 
than a systematic shift. It is interesting to note that in the RZA res- 
imulations, all the Virgo reproductions have a lower mass than the 
original BOX160 Virgo, just as the latter has a lower mass than 
th e observed Virgo clus ter. This may be connected to the findings 
of lCourtois et al.M2012h . who analysed the Cosmicflows-1 distance 
catalogue with the Wiener filter and found that the Virgo cluster is 
not dominating the peculiar velocity field as much as expected. It 
may therefore be harder to recover accurately in a constrained sim- 
ulation. 

Another possibility to estimate the robustness of the recon- 
struction is to compute an average of the density and velocity fields 
over the different evolved realizations. This way, structures com- 
ing from the random component will tend to average out and be 
suppressed, while structures appearing consistently in every real- 
ization will be enhanced. Figure [6] shows the same slice as in Fig- 
ure [5] but averaged over all six realizations A - F for the C30_10_I 
(without RZA), C30_10_II (with RZA), and E60_10_II (with RZA, 
larger mock). The position of the original BOX 160 Virgo cluster 
is marked with a blue cross. The resimulations done without RZA 
(top right) show some overdensity in this region, and a tendency 
of a flow towards it, but the overdensity is not high enough to 
create a cluster of mass comparable to Virgo in any of the non- 
RZA realizations. The density peak is much more pronounced in 
the simulations with RZA (bottom panels), which all have a mas- 
sive object near to that location. We also see that other structures 
with less mass, such as the overdense region below Virgo around 
X = 90 Y = 40, are not present in the averaged fields, which means 
that they lie in a mass range that is too low to be recovered by the 
reconstruction of initial conditions for any of the methods/mock 
catalogues. 



3.4 Positions and masses of other clusters 

To obtain a more precise estimate about the mass scale on which 
objects can be consistently recovered in the constrained resimula- 
tions using radial peculiar velocity constraints, we search for other 
massive clusters in the data zone. The largest overdense region 
within the data zone is dominated by the Centaurus and Hydra 
clusters, with masses of 6.06 x \0 u M Q /h and 5.20 x 10 14 M o //*, 
respectively. These two clusters are separated by a distance of 10 
Mpc/h. This overdense region is robustly recovered in all simula- 
tions including the ones from Method I not using RZA. It therefore 
lies above the minimum mass scale that can be recovered with- 
out RZA. Figure |7] shows the corresponding slice of the density 
and velocity fields that contains the Hydra and Centaurus clusters, 
averaged over the different realization^ All resimulations show a 



3 Note that in the observational data the Hydra/Centaurus clusters and the 
Virgo cluster all l ie within the supergalactic plane around SGZ = (see e.g. 



Courto is et al.l2012l) . On the other hand, while BOX 160 uses the same coor- 
dinate system, there the Hydra/Centaurus clusters and the Virgo cluster are 
not located at the same Z (=in the same X/Y plane), but actually in planes 
about 10-15 Mpc/h apart from each other in the Z dimension. This is an- 
other example of the systematic shifts that occur in constrained simulations 
from peculiar velocity data if one does not use Lagrangian reconstruction 
such as RZA. 
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Simulation 


Position X,Y,Z 


pec. vel. v x ,v y , v z 


Mass Position error 


Mass relative to 




[Mpc//z] 


[km/s] 


[10 u M Q /h] [Mpc/h] 


BOX160 Virgo 



BOX160 Virgo 82.0,74.4,91.8 +22, +142, -473 3.25 



C30_10_II_A 


84.8,78.5,92.6 


-23,-61,-391 


1.66 


5.0 


51 % 


C30_10_II_B 


80.5,76.8,96.7 


-22, +6, -912 


1.98 


5.7 


61 % 


C30_10_II_C 


80.3,72.6,94.4 


+26, +135, -403 


3.24 


3.6 


100% 


C30_10_II_D 


83.6,73.8,92.1 


+ 101, +231, -326 


1.21 


1.7 


37 % 


C30_10_II_E 


80.2,76.5,94.3 


+72, +247, -755 


2.33 


3.7 


71 % 


C30_10_II_F 


87.1,78.0, 92.1 


-273, +75, -527 


2.83 


6.3 


87 % 


C30_10_II mean 


82.8,76.0,93.7 


-19, +106, -552 


2.21 


2.6 


68% 


standard deviation 


2.6,2.1,1.6 


122,112,212 


0.68 


4.6 


21% 



E60_10_II_A 


82.8,78.5,90.5 


-113,-64, -507 


1.87 


4.4 


58% 


E60_10_II_B 


78.6, 74.7,93.9 


-200, +149, -682 


2.73 


4.1 


84% 


E60_10_II_C 


77.1,72.8,91.8 


+29, +135, -374 


1.21 


5.2 


37 % 


E60_10_II_D 


82.0, 76.0, 89.5 


-142, +201,-576 


3.03 


2.8 


93 % 


E60_10_II_E 


78.5,75.3,92.2 


+85, +289, -706 


1.31 


3.7 


40% 


E60_10_II_F 


84.0, 76.3,88.8 


+21, +259, -486 


2.80 


4.1 


86% 


E60_10_II mean 


80.5,75.6,91.1 


-53, +161, -555 


2.16 


2.1 


66% 


standard deviation 


2.5,1.7,1.7 


103,115,115 


0.73 


4.1 


22% 



Table 1. Virgo candidates found in the RZA resimulations of BOX 160 at z = 0. The first line corresponds to the original "Virgo" object in the BOX 160. The 
following lines list the haloes found in the AHF catalogues of the respective resmulations that are within 10 Mpc/h of the BOX 160 Virgo position and have a 
virial mass of at least 10 14 M o //l 



massive overdensity in this region. However in the Method I res- 
imulations there is a clear systematic shift in position. Further, in 
the averaged map one can see only one smeared out peak instead 
of two, which is additionally shifted towards positive X and nega- 
tive Y. In two of the six C30_10_I realizations (B and E), we can 
find only one massive cluster around the region where the Hydra- 
Centaurus pair should be. In these realizations, either the two over- 
dense regions merged during their evolution due to shifts in their 
position and/or displacement, or already in the reconstructed initial 
conditions the accuracy was not sufficient to robustly resolve them 
into two separate peaks. On the other hand, in the RZA resimula- 
tions, both clusters are resolved robustly and show up as separate 
peaks that are very close to their intended position; in every realiza- 
tion using Method II, we can find appropriate objects for Hydra and 
Centaurus at z = within a few Mpc/h of their original positions, 
just like for the Virgo cluster. Like in the case of Virgo, the masses 
of the resimulated Hydra/Centaurus clusters at z = show both a 
scatter and a systematic deviation from the original masses of the 
original BOX160 clusters of 6.06x \0 u M Q /h and 5.20x 10 14 M Q /h, 
respectively. The average mass of the resimulated Centaurus is at 
90% of the original mass in BOX160, with a scatter within a factor 
of two. Interestingly, the resimulated Hydra cluster is more mas- 
sive than it should be, at 173% of the original mass on average. In 
three out of the six C30_10_II realizations it is the more massive 
of the pair, although it should be the less massive, and in four out 
of six it even has a mass slightly above 10 15 M Q /h. As it is in the 
case of the Virgo cluster, this systematic error in the cluster masses 
does not improve if one goes from the C30_10_II realizations to the 
E60_10_II realizations. 

If we add more constraints and increase the data volume, 
the mass scale that can be reproduced with the RZA method in- 
creases noticeably. For the C30_10_II simulations, we can always 



find a Virgo cluster, but already the next massive cluster (the 
fourth-massive halo within 30 Mpc/h) with 2.41 x \0 u M Q /h can- 
not be unambiguously identified in one of the six realizations, the 
C30_10_II_E. Namely, within a distance of 10 Mpc/h and a mass 
within a factor of 3 there is no matching object in this realiza- 
tion. The next-massive clusters (ranked 5 to 7 by mass within 30 
Mpc/h) cannot be reliably found anymore. On the other hand, if 
we go to the E60_10_II simulations, we find resimulated counter- 
parts for all seven most massive clusters (the seventh having a mass 
of 0.96 X \0 u M Q /h) in all six realizations, all within a factor of 
3 in mass and within 7 Mpc/h in distance. The seventh most mas- 
sive halo with 0.96 x 10 14 M Q /h appears as a clear peak in the av- 
eraged map of E60_10_II realizations (Figure [7] bottom right, la- 
belled "cZ") close to its original position. For even lighter objects, 
it is not possible anymore to find unambiguous counterparts in the 
E60_10_II realizations. This is aggravated by the fact that below 
a certain mass there is a quickly increasing probability to find a 
seemingly matching, but actually randomly created object to ap- 
pear around the right position. 

The origin of the systematic errors on the reproduced ob- 
jects' masses can be understood from the non-linearity of the struc- 
ture formation process. The typical mass scatter within a factor 
of two that we observe here is consistent with the findings of 
iLudlow & Porcianil fooill) , who studied the relation between viri- 
alised haloes and their protohalo peaks in the initial conditions. 
They found that a protohalo peak on some fixed scale determines 
the mass of the resulting halo only within a factor of two. This 
explains the cluster mass discrepancies in our constrained resim- 
ulations. The exact virial mass of an object at redshift z = is a 
product of various non-linear structure formation processes, such 
as at what rate it can accrete mass from the surrounding structure 
and how efficiently it is fed from outside by connected filaments. 
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Figure 6. Density and velocity fields in a 10 Mpc/h thick slice at 87 < Z < 97 Mpc/h for the original BOX160 simulation (top left) and the average of all six 
constrained realizations without RZA (top right), with RZA (bottom left), and with RZA from the E60.10 mock using a larger data volume (bottom right). 
The density fields were smoothed with a Gaussian of radius 2.5 Mpc/h. The original centre of the BOX160's Virgo cluster is marked with a blue cross in each 
map. The thick contour line marks the cosmic mean density; solid contour lines are drawn for overdensities and dotted contour lines for underdensities. 



Such details of the structure around clusters cannot be recovered 
from reconstructed linear initial conditions. In the mass function 
of the 30 Mpc/h data zone (Figure |4j we see that BOX 160 has a 
very specific distribution: two objects (Hydra and Centaurus) in the 
highest-mass bin, one object (Virgo) in the next bin, and no objects 
in the bin after that. The average mass function of the constrained 
realizations does not follow this peculiar shape, but instead tends 
towards a distribution that is statistically more likely to occur. 



3.5 Filaments and voids 

Apart from massive dark matter haloes, another characteristic of 
the large-scale distribution of matter is the presence of filaments 
and voids. An example in the selected BOX160 data zone is the 
filament parallel to the F-axis, which goes upward from Virgo to- 
wards the BOX160's Coma cluster (outside of the map); it is not re- 
covered consistently in the RZA resimulations. We saw already in 
the RZA reconstruction of the displacement field (Paper I), that the 
RZA reconstruction struggles to recover such filamentary features 
accurately. This could be due to the non-linear structure formation 
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Figure 7. Same as Figure|6] but for a different slice at 69 < Z < 79 Mpc/h. In this slice, the BOX160 contains three massive clusters: Centaurus (Cen) with 
virial mass 6.06 x 10 14 M Q /h, Hydra (Hya) with 5.20 x 10 14 M o //z, and "Cluster Z" (cZ) with 0.96 x 10 14 M Q /h. The original centres of these three objects in 
BOX 160 are marked with blue crosses in each map. 



in such regions: in general, structure formation first proceeds along 
one dimension, forming sheets, followed by a collapse in the next 
dimension, forming fil aments, and finally the mater ial in these fila- 
ments falls into haloes dShandarin & Zeldovichl 1989b . The peculiar 
velocity field at z = only retains information about this last stage. 
In the particular case of the filament that connects the Virgo cluster 
to the Coma cluster in BOX 160, the input data contains only the 
infall of objects along the filament toward Virgo, but not the initial 
velocity distribution that led to the formation of that filament, and 
therefore our reconstruction procedure (WF + RZA) cannot recover 
it. In the case of the C30_10_II_A resimulation, this filament is in- 
stead replaced by another filament that is imposed by the random 



modes and aligned differently. This behaviour is repeatedly seen in 
the resimulations at other locations as well. At the other end of the 
filament there is a flow towards the Coma cluster seen at the up- 
per edge of the slice in Figure [5] (around X = 80; Y > 100). This 
flow is missed by all resimulations that use the C30_10 mock, since 
the BOX160 Coma cluster is too far away from the data zone. On 
the other hand, in the E60_10_II resimulations, this flow is recov- 
ered accurately because of the enlarged data zone, but this is not 
sufficient to also faithfully reproduce the alignment of the filament. 

It is not entirely clear at this point why the alignments of fil- 
aments are not always correctly reproduced in constrained simula- 
tions. Besides the inherent non-linearity of filament formation and 
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therefore the inability to extrapolate this process back in time with 
RZA another possible reason could be the fact that the input data 
consist of only radial peculiar velocities and the accuracy of fila- 
ment reconstruction could depend on how the filament is aligned 
compared with the line of sight. 

We therefore cannot generally trust that the alignment of fil- 
amentary structure in the Local Universe can be sufficiently re- 
produced in constrained simulations, except if their formation is 
strongly constrained by massive objects inside the data zone, as is 
the case for the thick filament hosting the Virgo cluster. The pres- 
ence of the latter is generally reproduced in all resimulations with 
RZA, although again its exact alignment is unconstrained. 

We did not compare the alignment of voids in the resimula- 
tions in detail, but we see that the presence and alignment of the 
most prominent voids in the data zone is generally recovered well. 
The largest voids develop by expansion of underdensities in the 
initial conditions that are above the minimum scale of initial condi- 
tions reconstruction and are therefore sufficiently constrained. The 
reconstruction of the displacement field using RZA helps track the 
expansion of voids back in time and generate appropriate initial 
conditions for them. Since it is believed that in the Local Universe, 
the outward push caused by the L ocal Void has an important in- 
fluence on shaping the Local Flow dTullv et alJl2008h . it would be 
interesting to study in detail to what extent this behaviour can be 
seen in the constrained simulations. This will be the subject of fu- 
ture work. 

3.6 The re-simulated peculiar velocity field 

To conclude the analysis of our re- simulations, we want to mention 
the reproduction quality of the BOX 160 peculiar velocity field at 
z = 0. We see that within the datazone, the peculiar velocity field is 
recovered exceptionally well in the RZA resimulations, especially 
for the better E6CL10 mock. We already saw that in this case the 
deviation of the resimulations from the original BOX 160 peculiar 
velocity field is only on the level of 1/4 of the total standard de- 
viation of the field. Comparing in Figure [5] the peculiar velocity 
field of the original BOX160 (upper left) and the E60_10_II res- 
imulation (bottom right), the structure of these two velocity fields 
is practically identical. The flows towards the major overdense re- 
gions and outwards of the large voids are all present with a correct 
alignment and amplitude. We can now argue that the limited ability 
to recover correct masses for particular clusters or the correct align- 
ment of filaments plays a lesser role and does not substantially limit 
the overall meaningfulness of the constrained simulation, since the 
large-scale velocity field is recovered accurately by the constrained 
resimulations. If we imagine that a Local Group-like object would 
be placed at the centre of the E60_10_II resimulation in Figure |5] 
then it would experience practically the same large-scale flow and 
would be embedded in a very similar large-scale environment com- 
pared to the original BOX160, so that we can study its dynamics. 
This is precisely the main motivation of running constrained sim- 
ulations, where of course the reference universe is the actually ob- 
served Local Universe. 



4 SUMMARY AND CONCLUSION 

In this paper, we investigated cosmological simulations that are 
designed to reproduce the large-scale structure of the Local Uni- 
verse. Such simulations can be produced by constraining their ini- 
tial conditions with the constrained realizations (CR) algorithm of 



Hoffman & Ribak (1991), using radial peculiar velocity data as in- 
put. On top of the previously established CR method we have added 
the RZA reconstruction, which is a novel Lagrangian reconstruc- 
tion scheme designed for peculiar velocity data. This new method 
allows to produce significantly more accurate constraints for the 
initial conditions. To quantify the accuracy of the method, we per- 
form tests using mock data extracted from a previous simulation 
serving as the "reference universe". We use different mock cata- 
logues with realistic observational errors that mimic the real data 
used for constrained simulations. After a reconstruction of the ini- 
tial conditions from the mock data, we evolve them forward again 
until z = in a series of re- simulations and compare their outcome 
to the original reference simulation at z = 0. We do so for both the 
previous method (without RZA) and our new method (with RZA). 

We find that with the previous method of generating con- 
strained initial conditions, if we use a sparse mock limited to 30 
Mpc/h and do not employ a Lagrangian reconstruction scheme (as 
in previous constrained simulations within the CLUES project), 
the threshold for robustly recovering structures is around ^ 5 x 
10 14 M o //*, and the accuracy on their positions is typically at the 
scale of 10 Mpc/h. Adding the RZA reconstruction, we can lower 
this threshold to under 3 x 10 14 M o //z, enough to robustly recover 
an object in our test simulation that is similar to the Virgo cluster. 
If we use the more complete E60_10 dataset, which covers a larger 
volume and approximates the accuracy of upcoming new datasets, 
the quality of reconstruction without using RZA does not increase 
much; but if we use RZA, we achieve a significantly better resolu- 
tion that allows us to robustly recover structure down to mass scales 
of about 1 x 10 14 M Q /h. Additionally, the RZA method reduces the 
errors on the positions where the resimulated structures appear to a 
scatter within about 5 Mpc/h around the true position. The average 
position of the different realization is even only 2 Mpc/h from the 
true position, meaning that there is practically no systematic shift 
of structures; this is a significant improvement over previous con- 
strained simulations from peculiar velocities, which featured shifts 
of the order of 10 Mpc/h and more due to a failure to account for 
the cosmological displacement field. 

We find that even with RZA and good input data, the method 
generally cannot reliably recover structures on mass scales below 
» 1 x \0 u M Q /h and also struggles to reproduce other large-scale 
features like the alignment of filaments. On the other hand, the pe- 
culiar velocity field at z = is reproduced exceptionally well by 
the re- simulations. This may be due to the fact that these velocities 
were used to constrain the realizations in the first place, and that 
they are less susceptible to non-linear effects than the cosmic den- 
sity distribution. The ability to reproduce the peculiar velocity field 
makes constrained simulations constructed with the RZA method 
an ideal laboratory to study velocity flows in the Local Universe 
when applied to actual observational data. We expect that with real 
data, the same degree of accuracy can be obtained with our tech- 
nique as we showed in the test presented here. Such an application 
of the method to the newest observational peculiar velocity data 
is a main focus of our future work. We expect that these simula- 
tions will present a significant methodical improvement over the 
currently available CLUES simulations and will provide a very use- 
ful framework for studying the dynamics of the Local Universe. 
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